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Abstract 


Statistical  resampling  plans  can  be  thought  of  as  procedures  for  "recycling  data."  They 
commonly  involve  computer-intensive  methods  for  the  interrogation  of  data  through  repeated 
inspection.  The  procedures  detailed  in  this  document  free  the  data  analyst  from  two  major 
restrictions  underlying  much  of  the  development  in  mathematical  statistics:  the  normal  or 
Gaussian  assumption,  and  the  necessity  to  focus  on  statistics  that  are  in  closed  form  and 
mathematically  tractable. 
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1.  The  Bootstrap 


Let  x  =  (xi,  X2,  •  •  • ,  xn)  be  a  random  sample  from  a  population  with  an  unspecified  distribution 
F  and  let  s(x)  be  an  arbitrary  statistic.  If  the  statistic  estimates  a  parameter  9  of  the  distribution 
F,  then  we  will  write  6  =  s(x),  where  the  hat  notation  denotes  an  estimate.  To  have  value  as  an 
estimator,  a  statistic  needs  to  estimate  accurately;  common  measures  of  accuracy  are  standard  error 
and  bias.  Almost  without  exception  (the  sample  mean  being  a  notable  exception)  an  estimate  of 
the  standard  error  of  a  statistic  cannot  be  determined  from  a  single  random  sample.  The  bootstrap 
procedure  attempts  to  overcome  this  impediment  through  intensive  computation.  Simply  put,  the 
bootstrap  procedure  attempts  to  provide  an  estimate  of  standard  error  of  an  arbitrary  statistic  9 
given  only  a  single  random  sample.  With  this  as  its  goal,  the  name  “bootstrap”  is  not  a  misnomer. 

A  fundamental  notion  (or  artifice)  is  that  of  a  bootstrap  sample.  A  bootstrap  sample  x*  = 
{xi,X2,  ■  ■ . ,  x*)  is  a  random  sample  of  size  n  drawn  with  replacement  from  the  class  of  n  objects 
making  up  the  random  sample  x  =  (xi,  ®2,  •  •  • ,  xn),  and  so  we  might  have  xf  =  xig,  x£  =  X7,  x£  = 
X3, . . . ,  x*  =  X5,  with  some  of  the  xi,  X2, . . . ,  xn  values  appearing  several  times  in  the  bootstrap 
sample  and,  others,  not  at  all. 

Why  would  anyone  want  to  engage  in  such  an  exercise? 

Kiefer  and  Wolfowitz  (1956)  established  for  a  univariate  random  sample,  x  =  (xi,X2, . . .  ,  xn), 
that  the  corresponding  empirical  distribution  F  is  the  nonparametric  maximum  likelihood  estimate 
of  the  unknown  distribution  F.  This  result  supports  the  construction  of  bootstrap  samples  in  one 
dimension,  but  extension  to  a  multivariate  situation  is  not  direct,  and  an  alternative  rationale  for 
bootstrap  sampling  is  needed. 

Toward  this  end,  consider  the  so-called  Dirac-comb  density  estimator  associated  with  a  one¬ 
dimensional  data  set  x  =  (xi,  X2, . . . ,  xn).  The  Dirac-comb  density  estimator  is  given  by 

•^(x)  =  ~£^x-x*)’  (x) 

where  the  expression  6(-)  is  a  Dirac  6  (or  unit  impulse)  function,  which  can  in  turn  be  expressed  as 

(2) 
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Expression  (1)  is  a  function  that  is  zero  everywhere  except  at  the  data  points.*  At  each  data 
point  the  density  is  unbounded;  the  assignment  of  probability  mass  1/n  at  that  location  is  merely 
suggestive  of  the  limiting  process.  As  a  nonparametric  density  estimator,  fs(x )  is  unappealing.  It 
dictates  that  any  subsequent  experiment  can  only  provide  the  data  points  already  observed,  and 
moreover,  those  points  must  have  equal  probability  of  occurrence. 

The  representation  fs(x)  is  not  without  value  in  some  situations;  for  example,  the  mean  of  fs(x) 
is 


roc  A  i  ” 

H6  =  /  xfs(x)dx  =  ~2^Xi  =  x 
J-c6  n  i=1 


(3) 


and 


cr|  =  /  (x-x)2fs(x)dx=-^2(xi-x)2  = 

J-oo  n  i=1 


(4) 


The  Dirac-comb  density  estimator  may  be  easily  extended  to  higher  dimensions.  For  a  sample 
X  =  (Xi,  X2, . . . ,  X„)  of  p-component  random  vectors,  fs(x)  becomes 


Mx)  =  ^s(x-xt)t 
n^i 


where 


(5) 


(6) 


and  where  Xj  is  the  j th  component  of  an  observation  X  =  (xi,  £2, . . . ,  xp). 


*The  unit  impulse  function  (2)  is  more  commonly  expressed  as 

S(x  -  Xi)  =  (  °*  X  *  Xi 

v  '  \  00,  X  —  Xi 


and 


/: 


S(x  —  Xi)dx  =  1. 
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Since  in  practice  bootstrap  sampling  is  carried  out  in  both  univariate  and  multivariate  situations 
without  regard  to  dimension,  viewing  the  bootstrap  procedure  as  based  on  a  Dirac-comb  estimator 
results  in  a  more  unified  theory. 

Through  a  straightforward  Monte  Carlo  procedure  in  which  many  samples  of  size  n  are  drawn 
with  replacement  from  F — bootstrap  samples — an  approximation  to  the  sampling  distribution  of  9 
is  produced.  In  order  for  the  statistical  accuracy  of  an  estimate  to  be  determined,  properties  of  its 
sampling  distribution  must  be  established.  For  each  bootstrap  sample  x^,  x£, . . . ,  x£  we  determine  a 
corresponding  9*  =  s(x*),  i  =  1,2 ,...,&.  This  empirical  approximation  to  the  sampling  distribution 
of  6  allows  an  assessment  of  standard  error  to  be  made. 

Some  attention  to  notation  is  in  order.  If  the  standard  error  of  9  is  denoted  as  ctf(9),  the 
bootstrap  estimate  of  standard  error  should  be  written  crp(9*),  where  the  star  (*)  notation  continues 
to  indicate  an  estimate  of  9  corresponding  to  a  bootstrap  sample  from  the  empirical  distribution 
F.  The  number  of  distinct  bootstrap  samples  from  a  data  set  of  n  distinct  values  can  be  shown  to 
be  m  =  C^n_1.  For  a  modest  sample  size,  say  n  =  10,  this  means  that  92,378  bootstrap  samples 
must  be  evaluated,  each  with  a  probability  of  occurrence  pj  given  by  a  multinomial  distribution,  in 
order  to  determine  ap{9*).  The  computational  burden  imposes  practical  limitations  on  the  sample 
size  n,  but  formally 


m 

»>(«*>  =  >[»(*,*) -»0]2}1/J,  (7) 

i= i 

where  s(-)  =  Pjsfej)>  and  x^,  x^, . . . ,  x*n  are  all  distinct,  is  the  appropriate  expression. 

Of  course  <Jp{9*)  is  in  general  not  equal  to  ap{9)  (the  value  we  really  seek),  but  it  is  the 
bootstrap  approximation.  Efron  and  Tibshirani  (1993)  refer  to  it  as  the  ideal  bootstrap  estimate 
of  standard  error.  As  suggested,  the  ideal  may  not  be  easily  attained  because  of  the  attendant 
combinatorial  explosion,  but  an  approximation  to  the  ideal  may  be  determined  via  the  following 
algorithm. 

The  bootstrap  algorithm. 

1.  Obtain  b  independent  bootstrap  samples  xj,  Xj, . . . ,  x£  by  sampling  with  replacement  from 
a  random  sample  x  =  (xi,  X2,  •  •  . ,  xn). 
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2.  Evaluate  6*  =  s(x*),  i  =  1, 2, . . . ,  b. 


3.  Calculate 


est  cf(h  =  {£[«?  -  «*(-)]V(l>  -  1)>1/2, 


i—l 


where  §*(-)  =  Ei=J*/b. 

Notice  in  step  3  of  the  algorithm  use  of  the  notation  est  crp{0*)\  an  alternate  expression,  <jp{d*), 
could  have  been  used. 

An  important  asymptotic  result,  stated  here  without  proof,  is  that  in  the  limit  as  b  — »  oo,  the 
bootstrap  estimate  of  standard  error  converges  in  probability  to  the  ideal  bootstrap  estimate: 


lim  estcrp(6*)  =  c p(0*).  (8) 

6— KX> 

This  result  provides  relief  from  the  tedius,  and  often  impractical,  attempt  at  enumeration  of  the 
distinct  bootstrap  configurations  needed  for  evaluation  of  equation  (7). 

Before  proceeding  further,  an  example  is  in  order. 

Example  1.  Measures  of  dispersion  for  bivariate  data  are  often  relevant  in  applied  problems, 
including  ballistic  studies,  where  the  ordered  pairs  (xi,j/i)  correspond  to  impact  locations  on  a 
target.  Measures  of  bivariate  dispersion  have  been  extensively  studied  (e.g.,  Grubbs  [1964]),  in¬ 
cluding  an  intriguing  statistic  known  as  the  extreme  spread  or  bivariate  range.  For  a  sample 
(xi,  j/i),  (x2, 2/2 ) ,  -  -  - ,  {xni  Vn)  the  extreme  spread  may  be  expressed  as 

R-=max.[(xi-Xj)2  +  (yi-yj)2]1/2  i,j  =  l,2,...,n.  (9) 

The  extreme  spread  is  simply  the  maximum  of  the  distances  between  all  CV;  pairs  of  points  in 
the  sample. 

This  statistic  is  intuitively  appealing  and  simple  to  evaluate,  and  yet  knowledge  of  its  distri¬ 
bution  remains  largely  confined  to  Monte  Carlo  studies.  This  holds  true  even  under  a  normal 
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assumption  for  the  variates  (rri,yj),  except  for  the  trivial  case,  n=2,  in  which  R  can  be  shown  to 
follow  a  chi  distribution. 

In  Table  1,  impact  coordinates  for  10  rounds  of  small  arms  ammunition  are  listed;  a  scatterplot 
of  the  data  appears  in  Figure  1. 


Table  1.  Impact  coordinates  for  10  rounds  of  small  arms  ammunition  fired  against  a 
vertical  target  (measurements  are  in  inches) 
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Figure  1.  Location  of  impact  for  10  rounds  of  small  arms  ammunition. 

Estimates  of  standard  error  for  the  extreme  spread  from  the  impact  data  in  Table  1  were 
computed,  following  the  bootstrap  algorithm.  Those  estimates  of  standard  error  for  values  of 
b=50(50)300  appear  in  Table  2. 

A  histogram  of  b  =  300  bootstrap  replications  is  shown  in  Figure  2.  Although  Cf0  =  45  distinct 
values  of  extreme  spread  are  possible,  after  300  replications,  only  14  have  been  realized.  This  is 
not  unexpected;  the  distribution  will  be  negatively  skewed  under  the  influence  of  the  max  operator 
appearing  in  statistic  (9) . 
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Table  2.  Bootstrap  estimates  of  standard  error  (300  bootstrap  replications  gave  the 
tabled  values  of  estcrp(Q*)) 


b:  50  100  150  200  250  300 

est  <TP(0*)  5.74  5.72  5.84  5.97  6.07  6.06 


8  10  12  14  16  18  20  22  24  26  28  30  32 


Figure  2.  A  histogram  of  b=300  bootstrap  replications  of  the  extreme  spread  for  the 
impact  data  in  Table  1. 

Guidance  regarding  the  number  of  bootstrap  replications  necessary  to  estimate  a  standard 
error  with  an  acceptable  degree  of  precision  is  limited  mostly  to  empirical  evidence.  This  evidence 
suggests  that  a  few  hundred  bootstrap  replications  is  usually  adequate,  although,  on  occasion,  a 
smaller  number  may  suffice.  A  definitive  answer  is  strongly  case-dependent,  and  involves  the  rate 
of  convergence  in  expression  (8). 


2.  The  Parametric  Bootstrap 

There  is  no  inherent  requirement  that  the  bootstrap  approach  be  restricted  to  a  nonparametric 
setting.  If  there  is  support  for  a  specific  form  for  the  distribution  F,  then  F  can  be  estimated 
parametrically.  We  will  adopt  the  notation  Ff  to  indicate  that  a  specific  functional  form,  /(•),  has 
been  assumed  for  the  random  variable.  Suppose  for  the  impact  data  in  section  1  the  sample  z  = 
{(xi,  yi),  (x2, 2/2)) . . . ,  (xn)  yn )}  is  assumed  to  have  come  from  a  normal  population  with  parameters 
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and  £.  Estimates  of  the  mean  vector  and  covariance  matrix  will  be  given  by  (x,  y)f  and 


£(xi  -  x)2  E(xi  -  x)(yi  -  y) 
£(xi  -  x)(yi  -  y)  E (y*  -  y)2 


From  the  impact  data  in  Table  1,  we  obtain  a  mean  vector  jj,  =  (10.45, 14.10)'  and  covariance 


matrix  E  = 


27.74  -15.27 

-15.27  67.76 


The  bootstrap  procedure,  after  parameterizing  F  (the  empirical 


distribution  function  F  is  replaced  by  a  normal  distribution  Fj  =  N(jl,  E))  reduces  to  an  ordinary 
Monte  Carlo  procedure.  Instead  of  sampling  with  replacement  from  z  to  generate  bootstrap  samples 
z\,  Zj, . . . ,  z £,  Monte  Carlo  samples  are  drawn  directly  from  Ff  to  produce  z*,  i  =  1, 2, . . . ,  b.  The 
histogram  of  R*  values  for  b=300  is  shown  in  Figure  3.  The  estimate  of  standard  error  from  these 
data,  est  ap^(§*)  =  6.71 — compared  to  6.06  for  the  nonparametric  model. 


Taylor  and  Grubbs  (1975)  carried  out  extensive  Monte  Carlo  simulation  studies  of  the  extreme 
spread  distribution  for  small  sample  sizes,  under  a  standard  normal  assumption  for  the  underlying 
distribution  F,  and  concluded  that  the  distribution  is  adequately  described  by  a  chi  variate. 


Figure  3.  Histogram  of  R*  values  for  b=300  replications  under  a  normal  assumption: 

F  =  N(n,  E). 
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3.  Bias  Estimation 


The  bootstrap  procedure  described  in  section  1  is  not  limited  to  estimation  of  the  standard 
error  of  a  statistic  9  =  s(x);  it  can  be  used  with  only  minor  modification  to  provide  an  estimate  of 
bias.  The  bias  of  an  estimator  9  is  the  difference  between  the  expectation  of  9  and  the  parameter 
9  that  it  purports  to  estimate;  formally, 


bias  =  Ef(9)  —  9.  (10) 

An  unbiased  estimator  (i.e.,  a  statistic  9  satisfying  bias  =  Ef(9)  —  9  =  0)  is  desirable  inasmuch 
as  it  provides  a  measure  of  assurance  that  the  estimator  is  not  producing  values  that  are  either 
systematically  too  small  or  too  large. 

The  fundamental  problem  remains  unchanged.  Only  a  single  random  sample  x  =  (rci,  X2,  ■  ■  ■ ,  xn) 
from  an  unspecified  distribution  F  is  available.  Following  the  same  rationale  used  in  assessing  the 
standard  error,  the  ideal  bootstrap  estimate  of  bias  becomes 


bias  =  Ep(9*)  -  9 ,  (11) 

where,  as  before,  F  is  the  empirical  distribution  of  F,  9*  is  the  bootstrap  estimate  of  9,  and  9  is  the 
estimate  of  9  based  on  the  random  sample  x.  The  same  computational  burden  attendant  to  the 
evaluation  of  crp(9*)  in  section  1  is  now  imposed  on  the  determination  of  the  expectation  Ep{9*), 
and  is  addressed  in  much  the  same  manner.  Namely,  the  further  approximation 

_ _  i  J> 

est  bias  =  —  9*  —  9  (12) 

°  i= 1 

is  invoked.  To  calculate  est  bias,  it  is  only  necessary  to  incorporate  expression  (12)  at  the  final  step 
of  the  bootstrap  algorithm,  and,  if  desired,  est  bias  and  estop  can  be  determined  simultaneously. 

The  bias  of  R  based  on  b  =  300  bootstrap  replications  from  the  impact  data  in  Table  1  was 
determined  to  be  est  bias  =  —5.55. 

Should  we  have  much  confidence  in  this  estimate  of  bias?...  No. 
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The  sign  of  the  resultant  estimate,  if  not  the  magnitude,  was  never  in  doubt.  This  is  an  example 
in  which  the  bootstrap  fails. 

To  see  why  this  is  so,  consider  the  estimate  of  bias  (12).  In  that  expression  we  have  9  =  R(x),  the 
extreme  spread  of  the  basic  random  sample  x.  Clearly,  no  bootstrap  sample  can  lead  to  a  9*  larger 
than  9.  It  will,  however,  produce  estimates  6*  less  than  9  whenever  a  pair  of  impact  coordinates 
maximizing  expression  (9)  fail  to  appear  in  the  bootstrap  sample  x*.  Thus  the  bootstrap  estimate 
of  bias  will  always  be  negative,  except  in  the  rare  event  when  a  maximizing  pair  appears  in  every 
bootstrap  sample. 

The  situation  is  closely  related  to  the  well-known  problem  of  estimating  the  range  of  a  uni¬ 
form  variate  C/[0, 9}  from  a  random  sample  x  =  (xi,  X2, . . . ,  xn),  in  which  the  maximum  likelihood 
estimate  is  9  =  X[n],  the  largest  order  statistic. 

A  reality  of  the  bootstrap  is  that  the  procedure  is  held  captive  in  the  convex  hull  of  the  original 
sample  x  =  (xi,X2, . .  • , xn).  A  modification  known  as  the  smoothed  bootstrap  has  been  advanced, 
in  which  a  scaled  random  variate  with  small  dispersion  is  added  to  each  observation  in  the  original 
data  set,  but  the  result  is  essentially  a  perturbed  Dirac-comb  density  estimator. 

This  is  not  a  contrived  example.  Instances  in  which  Dirac-comb  density  estimator  based  tech¬ 
niques  (viz.,  the  bootstrap)  produce  unsatisfactory  results  are  not  difficult  to  find.  Consider  the 
situation  where  several  rounds  of  ammunition  impact  nearby,  but  never  directly  on,  an  aim  point.  If 
the  distribution  of  shots  is  assumed  circular  normal  with  mean  located  at  the  aim  point — the  most 
common  assumption  in  ballistic  analyses — then  a  neighborhood  of  the  aim  point  will  appear,  with 
high  probability,  to  be  a  particularly  safe  place  to  reside.  Anomalies  such  as  these  gave  impetus  to 
an  algorithm  advanced  by  Taylor  and  Thompson  (1982,  1986),  which  is  discussed  in  section  5. 

4.  Bootstrap  Confidence  Intervals 

Application  of  the  bootstrap  procedure  to  confidence  interval  construction  has  been  a  major 
research  focus,  and  a  substantial  number  of  procedures  toward  this  end  have  been  advanced.  In  this 
section  a  single  approach,  the  percentile  method,  for  constructing  approximate  confidence  intervals 
in  small  sample  nonparametric  situations  is  considered.  The  intent  is  not  to  suggest  that  this  is  in 
some  sense  a  “best”  method;  it  is  rather  a  widely  used  approach  whose  motivation  is  reasonably 
accessible. 
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Under  quite  general  conditions,  the  distribution  of  an  estimate  9  of  a  parameter  9  approaches  a 
normal  distribution  as  the  sample  size  x  =  ( xi,X2 ,  ■  . .  ,xn)  increases  without  bound.  This  asymp¬ 
totic  result  leads  to  the  approximation 


9  ~N{9,  a2) 

from  which  a  100(1  —  2a)  confidence  interval  for  6. 

[6  —  z (1-Q)  •  a,  6  —  z ^  •  a] 


(13) 


(14) 


commonly  expressed  as 


0±z(1-q)-ct  (15) 

is  derived.  This  is  sometimes  referred  to  as  the  standard  confidence  interval  for  0;  is  the 

100(1  —  £*)th  percentile  of  the  standard  normal  distribution. 

Now  consider  as  an  estimate  of  6,  the  bootstrap  estimate,  0*.  Appealing  to  the  same  argument, 
we  are  led  to  the  approximation 


0*~iV(0,a2),  (16) 

where  9  is  determined  from  x  =  (aq,  X2,  ■  ■  ■  ,xn)  and  a1  from  any  appropriate  estimation  procedure 
(e.g.,  the  bootstrap,  jackknife,  or  historical  evidence). 

Proceeding  formally, 


pr[2(«)  <  _  —  <  2;(1  “)] 

Pr[z{a)  •  a  <  0*  —  0  <  z^-Q)  ■  a) 
Pr[0  +  z(a)  -a  <9*  <9  +  z(1_a)  •  a] 


1  -2a. 


(17) 
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Notice  that  the  lower  and  upper  confidence  limits  for  6  appearing  in  (14)  coincide  with  the 
lOOath  and  100(1— a)th  percentiles  of  the  distribution  of  0*  (since  This  observation 

is  important  because,  when  the  distribution  of  0*  is  approximately  normal,  one  might  entertain  as 
an  approximate  1  —2a  central  confidence  interval  for  0  the  interval  [$*(“),  <?*(1-a)] — the  lOOath  and 
100(1  —  a)th  percentiles  of  0*. 

And  this  is  precisely  what  the  percentile  method  does. 

Further  investigation,  beyond  the  scope  of  this  exposition,  suggests  that  this  procedure,  or  a 
modification  of  this  procedure,  may  provide  acceptable  results  in  those  situations  for  which  there 
exists  a  transformation  to  normality  for  the  estimator  6*.  In  that  eventuality,  the  percentile  method 
exploits  this  circumstance  automatically,  and  without  the  necessity  of  formally  determining  the 
transformation. 

Consider  the  following  example. 

Example  2.  An  historic  approach  to  vulnerability  assessment  of  military  systems  commences  with 
the  overlay  of  an  mxn  rectangular  grid  as  illustrated  in  Figure  4.  For  each  cell  in  the  grid,  a  detailed 
computer  model  estimates  the  probability  of  incapacitation  of  the  system  that  would  result  from 
an  impact  in  that  location  under  prescribed  ballistic  conditions.  These  conditional  probabilities, 
commonly  referred  to  as  probability-of-kill,  will  be  denoted  by  ,  ?  =  1,2,...  ,mn.  The  cell-level 
P*.|/iis  are  combined  under  a  weighting  scheme  over  the  entire  grid  to  produce  a  global  estimate  of 
incapacitation, 


Figure  4.  Armored  vehicle  with  an  overlaid  rectangular  grid. 


Taylor  and  Bodt  (1984)  constructed  approximate  confidence  intervals  for  P*,,  using  the  percentile 
method,  for  several  choices  of  grid  size.  Coarse  grid  templates  were  chosen,  but  with  the  finest  level 
of  detail  retained,  to  produce  within-cell  variation.  A  set  of  Pk\hi  values  was  obtained,  one  value 
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from  each  grid  cell,  to  form  the  base  sample  from  which  bootstrap  samples  would  be  drawn.  Each 
confidence  interval,  [_p*(Q),  P*'1  °^],  was  based  on  b  =  1,000  bootstrap  estimates  of  Pk- 

This  approach  differs  from  the  standard  bootstrap  procedure  outlined  in  section  1  in  that  a 
stratified  sample,  rather  than  a  simple  random  sample,  was  collected.  The  modification  is  ap¬ 
propriate  for  this  particular  application,  and,  since  Pk  is  essentially  a  weighted  sum,  the  Central 
Limit  Theorem  supports  the  anticipation  of  an  approximate  normal  distribution  for  the  bootstrap 
estimates  P£. 

A  total  of  100  confidence  intervals  were  constructed  for  each  grid  size  by  returning  to  the 
coarsened  grid  100  times  to  establish  a  new  stratified  base  sample,  and  then  repeating  the  entire 
bootstrap  procedure.  This  exercise  is  again  outside  the  usual  bootstrap  venue;  the  intent  was 
to  see  if  one  hundred  (1  -  2a)-level  confidence  intervals  covered  the  parameter  Pk  approximately 
100(1  —  2a)  times.  Empirical  coverage  results  are  shown  in  Table  3  and  the  mean  interval  width  in 
Table  4. 

Table  3.  Number  of  bootstrap  confidence  intervals  covering  the  finest-level 

parameter  estimate,  Pk  (100  confidence  intervals  were  constructed  for  each 
table  entry) 


grid 

.99 

.95 

.90 

.80 

24x60 

100 

100 

100 

100 

16x40 

100 

100 

100 

99 

12x30 

100 

100 

100 

100 

8x20 

100 

99 

99 

95 

Table  4.  Mean  width  of  the  confidence  intervals  corresponding  to  the  entries  in 
Table  3 


grid 

.99 

.95 

.90 

.80 

24x60 

.056 

.042 

.035 

.027 

16x40 

.084 

.063 

.053 

.041 

12x30 

.112 

.084 

.071 

.055 

8x20 

.167 

.127 

.106 

.083 

Here  the  percentile  method  appears  overly  conservative.  Efron  and  Tibshirani  (1993)  suggest 
modifications  intended  to  increase  its  efficacy.  DiCiccio  and  Romano  (1988),  Hall  (1988),  and  Loh 
(1991)  provide  a  more  complete  discussion  of  confidence  interval  construction. 


12 


5.  A  Generalized  Bootstrap 


The  nonparametric  bootstrap  is  not  concerned  with  the  construction  of  a  nonparametric  density 
estimator  beyond  the  rudimentary  Dirac-comb:  its  focus  is  on  the  generation  of  values  that  appear 
to  conform  to  the  same  structure  as  that  underlying  the  authentic  data  set.  Given  the  ubiquity 
of  fast  computing,  consideration  of  resampling  schemes  based  on  estimators  more  complex  than 
the  Dirac-comb  appears  warranted,  and  that  is  the  focus  of  this  section.  A  resampling  plan  that 
is  useful  for  simulating  data  (SimDat)  from  an  unspecified  multivariate  distribution  is  developed. 
The  procedure  is  shown  to  behave  in  some  respects  like  a  normal  kernel  estimator;  i.e.,  the  density 
may  be  approximated  by 


h{X)  =  -YjK{x~Xu  Si), 

Tt  .  ^ 
i=l 


(18) 


where  K(')  is  a  normal  distribution  centered  at  zero  with  a  locally  estimated  covariance  matrix  E*. 
The  bootstrap  procedure  arises  as  a  limiting  case,  and,  in  that  sense,  SimDat  may  be  considered  a 
generalized  bootstrap. 

The  algorithm  operates  as  follows. 

Let  X  =  (Xi,  X2, . . . ,  Xn)  be  n  observations  of  a  p-component  random  vector  X  from  a  pop¬ 
ulation  with  an  unspecified  distribution  F.  Suppose  that  a  pseudo-sample  of  size  N,  where  n  and 
N  are  not  necessarily  equal,  is  to  be  generated.  The  data  are  rescaled  so  that  the  marginal  sample 
variances  coincide  and,  for  each  of  the  n  observations,  a  specified  number  (m-1)  of  nearest  neighbors 
is  determined.  The  nearest  neighbors  are  retained  in  an  nx(m-l)  array  for  use  as  a  reference  set. 

To  commence,  a  single  point  is  selected  at  random  from  the  data  set  X  and  its  m-1  nearest 
neighbors  recalled  from  the  reference  array.  The  mean  of  the  resultant  set  of  m  points, 

1  m 

*  =  (19) 

is  determined,  and  each  of  the  m  points  coded  about  X  to  obtain 

X[  =  Xi-X)  t  =  1,2,  (20) 
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This  calculation,  like  the  determination  of  the  m-1  nearest  neighbors  for  each  data  point,  will  be 
required  only  once.  The  X ■  values  along  with  X  will  be  retained  in  an  nx(m+l)  array. 

Next,  a  random  sample  «i,  «2>  •  •  •  > um  is  generated  from  the  uniform  distribution 


/3(m  —  1) 
m2  ’ 


1 

- h 

m 


3  (m  —  1)  \ 
m2  I 


and  the  linear  combination 


(21) 


X'  =  Y,uiX'i  (22) 

t=i 

formed.  Finally,  with  the  addition  of  X,  a  pseudo-data-point  X  is  produced: 


X  =  X'  +  X. 


(23) 


And  the  process  repeats. 

A  second  point  is  selected  at  random  from  X  and  its  m-1  nearest  neighbors  recalled;  continuing 
until  at  last  N  pseudo- values  have  been  determined. 

With  increasing  m  and  n,  the  procedure  behaves  very  much  like  the  normal  kernel  approach 
mentioned  earlier.  To  see  why  this  is  so,  consider  the  sampled  vector  X\  and  its  m-1  nearest 
neighbors: 


Xi  = 


xu 

X2  i 


|_  Xpi 


,  i  =  l,2,...,m. 


(24) 


For  the  moment,  consider  this  collection  of  m  points  as  coming  from  a  truncated  distribution 
with  mean  vector  ft  and  covariance  matrix  E.  For  the  random  sample  ui,U2,  ■ . . ,  um  from  the 
uniform  distribution  (21),  we  have 


.  1  T.  .  .  m  —  1 

E[ui)  —  — ,  Var\Ui)  —  2  > 


m 


m* 


(25) 
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and 


Cov(ui,  Uj)  =  0,  for  i  ^  j. 


(26) 


In  the  linear  combination 


m 

Z  =  (27) 

t= 1 

the  rth  component  of  the  vector  Z,  zr  =  u\xr\  +  U2Xr2  H - h  umxTm ,  satisfies 


—  Mr  j 


(28) 


Var{zT )  =  of  + 


m  —  1 
m 


(29) 


and 


Cov(zr,  zs)  =  ars  + 


m  —  1 


m 


Mr  Ms* 


(30) 


If  the  mean  vector  of  the  random  variable  X  was  0  =  (0, 0, ,  0)',  then  the  mean  vector 
and  covariance  matrix  of  Z  would  be  the  same  as  that  of  X:  E{zr)  =  0,  Var{zT)  =  of,  and 
Cov(zr,zs )  =  ars.  By  translation  to  the  local  sample  mean  of  the  nearest-neighbor  cloud — the 
motivation  behind  expression  (23) — we  do  not  achieve  this  result  exactly,  but  we  do  come  very 
close  to  generating  an  observation  from  the  truncated  distribution  that  models  the  points  in  the 
nearest-neighbor  cloud.  For  moderately  large  m,  the  Central  Limit  Theorem  becomes  dominant, 
and  SimDat  comes  close  to  sampling  from  n  normal  distributions  with  mean  and  covariance  matrices 
corresponding  to  those  of  the  n  m-nearest-neighbor  clouds. 

Sufficient  conditions  for  consistency  of  the  nonparametric  density  estimator  corresponding  to 
SimDat  is  problematic.  Mack  and  Rosenblatt  [10]  advance  a  formula  for  nearest-neighbor  nonpara¬ 
metric  density  estimators: 


m  =  Cn4/(P+ 4)f  (31) 


which  seems  promising,  were  it  not  for  the  fact  that  the  parameter  C  is  usually  not  available. 
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At  the  extremes,  a  choice  of  m=l  causes  SimDat  to  collapse  to  the  bootstrap;  m=n  produces  an 
estimator  that  roughly  samples  from  a  multivariate  normal  distribution  with  the  mean  vector  and 
covariance  matrix  computed  from  the  original  data  set. 


* 


* 

+ 

*  +  **  +  # + 

*  +  *  +  *,+ 

*  *  T  + 

♦  4.  *  " 

*  + 

■*  * 

**+  \  * 

* 

*  * 

+  * 

* 

1  1  1  1  *' 

*  +*  ^A***<\, 
*  *  3*  **  * 

* 

±  * 

V* 

*  +  + 

*  V  * 

* 

¥ 

the  actual  data 

m=5 

+  is  the  simulated  data 


Figure  5.  Samples  from  a  normal  mixture:  l/2Ni  +  I/3N2  +  I/6N3  with  m=5. 

Example  3.  In  Figure  5,  a  sample  of  85  values  from  a  mixture  of  3  normal  distributions, 
Ei),  i  =  1,2, 3,  with  mean  and  covariance  matrices 


Mi 


M2  = 


M 

,  £1  = 

^  1  1/2  ' 

i-2/ 

^-1/2 1  j 

,  £2= 

1 1/2  )• 

00 

(1/2 1 ) 

and 


(32) 


(33) 


M3 


3/2  ) 


£3  = 


V  1/10 


(34) 


along  with  a  pseudo-data-set  of  size  85  generated  by  SimDat,  with  m=5,  is  presented.  The  emulation 
of  the  data  appears  satisfactory,  at  least  visually.  In  Figure  6,  the  same  exercise,  but  with  111= 15, 
is  shown.  There,  effects  of  a  modest  oversmoothing  appear  to  be  present. 
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*  is  the  actual  data 


m=15 


+  is  the  simulated  data 

Figure  6.  Samples  from  a  normal  mixture:  l/2Ni  +  I/3N2  +  I/6N3  with  m=15. 


Example  4.  The  SimDat  algorithm  was  applied  to  the  impact  data  from  Table  l.t  Three  hun¬ 
dred  pseudo-data-sets  of  size  10  were  generated,  with  a  value  of  m=3,  and  their  extreme  spreads 
determined.  A  histogram  of  the  results  is  shown  in  Figure  7. 


Figure  7.  SimDat  algorithm  applied  to  the  impact  data,  with  m=3. 

Empirical  evidence  suggests  that  if  the  data  set  is  very  large,  say  of  size  1,000  or  greater,  good 
results  are  generally  obtained  with  m«.02n.  For  smaller  values  of  n,  m«.05n  often  work  well. 

t  A  version  of  SimDat  in  the  S  language  is  available  under  the  name  “gendat”  from  the  S  Library;  a  Fortran  version 
of  SimDat  entitled  “RNDAT”  is  available  from  IMSL. 


6.  Summary 


A  comparison  of  the  nonparametric  bootstrap,  parametric  bootstrap,  and  SimDat  procedure, 
all  applied  to  estimation  of  extreme  spread  from  the  impact  data  presented  in  Table  1,  is  shown 
in  Table  5.  This  application  proves  vexing  to  any  approach,  as  the  tabled  values  suggest.  The 
nonparametric  bootstrap  applied  to  the  impact  data  has  been  extensively  discussed.  It  is  worth 
remembering  that  the  extreme  spread  of  the  original  data  set,  9  =  30.21  (appearing  as  the  maximum 
value  for  the  nonparametric  bootstrap)  is  the  sole  authentic  estimate.  The  parametric  bootstrap 
has  the  smallest  estimate  of  bias,  a  direct  consequence  of  the  strong  distribution  assumption.  Notice 
the  maximum  value  of  54.64 — far  in  excess  of  30.21.  SimDat  has  the  smallest  estimated  standard 
error  and  the  largest  bias.  It  generated  some  extreme  spreads  larger  than  the  authentic  6 — notice 
the  maximum,  31.18 — but  most  pseudo-points  are  aggregated  in  the  regions  of  highest  population 
density  by  virtue  of  the  randomly  selected  center  points  that  commence  the  algorithm.  In  other 
words,  SimDat  is  providing  a  modest  amount  of  smoothing  under  control  of  the  parameter  m,  as 
it  was  intended  to  do. 

Table  5.  Descriptive  statistics  for  the  nonparametric(n)  and  parametric(p)  bootstrap 
and  the  SimDat  procedure.  Estimates  of  standard  error,  bias,  and  maximum 
and  minimum  values  of  the  extreme  spread  are  for  300  replications. 


method 

std  err  est 

bias  est 

min 

max 

bootstrap(n) 

6.06 

-5.55 

9.82 

30.21 

bootstrap(p) 

6.71 

-2.60 

13.18 

54.64 

SimDat 

4.71 

-11.56 

7.78 

31.18 

And  finally,  a  note  of  caution.  Both  the  nonparametric  bootstrap  and  SimDat  are  model-free 
techniques  for  examining  data.  The  nonparametric  techniques  use  the  power  of  the  computer 
to  circumvent  detailed  modeling,  and  such  an  approach  is  highly  useful  as  an  interpolative  device. 
When  the  dimensionality  of  a  data  set  becomes  high-say,  five  or  greater,  these  approaches  should  be 
used  with  caution,  since  a  number  of  widely  separated  modes,  with  data-scarce  regions  in  between, 
might  easily  be  obscured.  Dealing  with  such  data  sets  nonparametrically  and  away  from  the  modes 
becomes  more  of  an  extrapolation  problem,  and  blindly  using  smoothed  interpolation  routines  can 
be  a  precursor  to  faulty  inference. 
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